function G = monitor(x,y,grad)
%  function G = monitor(x,y,grad)
%   grad�Ǵ�СΪ��row,col������ԭͼ���Сһ��
%   x,y�ϴ֣�Ϊ��������
%   G��С��x,yһ�£����ݶ��ڸõ�Ԫ�ϵ�ƽ��ֵ����

% xmax = max(max(x));
% ymax = max(max(y));
[row,col] = size(x);
G = zeros(row,col);
row = row - 1;
col = col - 1;

% ���������gradΪÿ��Ԫ����montecarlo���,��Ҫ����֧�֣���
npt = 500;
for k = 1:npt
    xi = rand(row,col);
    eta = rand(row,col);
    xi = abs(xi - 1e-2);
    eta = abs(eta - 1e-2);
    a1 = (1-xi).*(1-eta);
    a2 = xi.*(1-eta);
    a3 = xi.*eta;
    a4 = (1-xi).*eta;
    z = x(1:row,1:col).*a1 + x(2:row+1,1:col).*a2 + ...
        x(2:row+1,2:col+1).*a3 + x(1:row,2:col+1).*a4;
    r = y(1:row,1:col).*a1 + y(2:row+1,1:col).*a2 + ...
        y(2:row+1,2:col+1).*a3 + y(1:row,2:col+1).*a4;
    idx1 = sub2ind(size(grad),ceil(z)  ,ceil(r)  );
    idx2 = sub2ind(size(grad),ceil(z)+1,ceil(r)  );
    idx3 = sub2ind(size(grad),ceil(z)+1,ceil(r)+1);
    idx4 = sub2ind(size(grad),ceil(z)  ,ceil(r)+1);
    xi = z - floor(z);
    eta = r - floor(r);
    a1 = (1-xi).*(1-eta);
    a2 = xi.*(1-eta);
    a3 = xi.*eta;
    a4 = (1-xi).*eta;
    val = grad(idx1).*a1 + grad(idx2).*a2 + grad(idx3).*a3 + grad(idx4).*a4;
    G(1:row,1:col) = G(1:row,1:col) + val;
end
G = G/npt;

G(:,col+1) = G(:,col);
G(row+1,:) = G(row,:);

% ���Ժ���1
%G = sin(pi*x/xmax).*sin(pi*y/ymax);

% ���Ժ���2
%G = 1.0 + 100*sin(pi*x/xmax).*sin(pi*y/ymax);

% ���ƺ���⻬
for k=1:3
    G(2:row,2:col) = 0.125*(4.0* G(2:row,2:col) + ...
         G(1:row-1,2:col) +  G(3:row+1,2:col) + ...
          G(2:row,1:col-1) +  G(2:row,3:col+1));
end

H = fspecial('gaussian',5,0.5);
%H = [0.011344 0.083820 0.011344; ...
%     0.083820 0.619347 0.083820; ...
%     0.011344 0.083820 0.011344];
G1 = conv2(G, H, 'same');
G1 = conv2(G1, H, 'same');

G0 = G/max(max(G));
G1 = G1/max(max(G1));

G = 1.0 + 1e7*G0 + 1e3*G1;
